/* normalgenerator.h
 *
 *	Copyright (C) 2008 CRIMERE
 *	Copyright (C) 2008 Jean-Marc Mercier
 *	
 *	This file is part of OTS (optimally transported schemes), an open-source library
 *	dedicated to scientific computing. http://code.google.com/p/optimally-transported-schemes/
 *
 *	CRIMERE makes no representations about the suitability of this
 *	software for any purpose. It is provided "as is" without express or
 *	implied warranty.
 *
 *  This program is free software: you can redistribute it and/or modify
 *  it under the terms of the GNU General Public License as published by
 *  the Free Software Foundation, either version 3 of the License, or
 *  (at your option) any later version.
 *
 *  This program is distributed in the hope that it will be useful,
 *  but WITHOUT ANY WARRANTY; without even the implied warranty of
 *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 *  GNU General Public License for more details.
 *
 *   You should have received a copy of the GNU General Public License
 *   along with this program.  If not, see <http://www.gnu.org/licenses/>.
 *
 *
 */
#include <src/algorithm/random/boost_sobol.h>
#include <src/math/distributions/distributions.h>
#include <src/miscellaneous/OS_basic.h>


#if !defined(_NormalGenerator_)
#define _NormalGenerator_
namespace boost {
namespace random {
/*! \brief Generator helper for the OS_inverse_cumulative_normal_distribution function.
	Given a uniform generator G, return a value corresponding to a normal distributed sampling :
	InverseCumulativeGenerator() return N^{-1}(G), where G++ is the next sampling of the uniform generator, and N 
	the gaussian distribution.
*/
template <class UniformGenerator_type = sobol_generator<32>, std::size_t Dim = 1>
class GaussianGenerator 
 {
	 OS_STATIC_CHECK(false); // to implement / to detail
 };

template <class UniformGenerator_type>
class GaussianGenerator<UniformGenerator_type,2> 
 {
 public:
	typedef GaussianGenerator type;
	typedef typename UniformGenerator_type::result_type result_type;
	GaussianGenerator(const result_type average, OS_double sigma, const UniformGenerator_type & u) : average_(average), Gaussian_(0.,sigma), UniformGenerator_(u)
	{};
	OS_double operator()() // next sampling
	{
		return Gaussian_(jmmath::norm2(UniformGenerator_()-average_));
	};

	template <class T, std::size_t Dim = T::Dim>
	struct fill_helper
	{
		fill_helper(type & c) : c_(c) {};
		void operator()(T & t) {
			typedef T::iterator iterator;
			typedef T::reference reference;
			iterator beg = t.begin();
			iterator end = t.end();
			for (;beg != end;++beg) {
				fill_helper<reference,Dim-1> test(c_);
				test(*beg);
			};
		};
		type & c_;
	};
	template <class T>
	struct fill_helper<T,1>
	{
		fill_helper(type & c) : c_(c) {};
		void operator()(T & t) {
			typedef T::iterator iterator;
			typedef T::value_type value_type;
			iterator beg = t.begin();
			iterator end = t.end();
			for (;beg != end;++beg) {
				*beg = c_();
			};
		};
		type & c_;
	};


	template <class T>
	void fill(T &t) {
		fill_helper<T>(*this)(t);
	};

 private :
	result_type average_;
	UniformGenerator_type UniformGenerator_;
	base_normal_distribution<OS_double> Gaussian_;
 };


template <class UniformGenerator_type>
class GaussianGenerator<UniformGenerator_type,1> : public Generator_facilities<GaussianGenerator<UniformGenerator_type,1>>
 {
 public:
	typedef OS_double result_type;
	GaussianGenerator(
		OS_double average = 0.0, 
		OS_double sigma   = 1.0, 
		const UniformGenerator_type & init = 0 ) :
	Gaussian_(average,sigma)
	{
		if (&init)
			UniformGenerator_ = init;
	};
	OS_double operator()() // next sampling
	{
		double test = UniformGenerator_();
		return Gaussian_(test);
	};
 private : 
	 UniformGenerator_type UniformGenerator_;
	 base_normal_distribution<OS_double> Gaussian_;
 };

} // namespace random
} // namespace boost


#endif
